www.gusucode.com > C++ 图像分割源代码 > C++ 图像分割源代码/gusucode/seg_source/Analysis.cpp

    //Download by http://www.NewXing.com
//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#include "./Analysis.h"
#include "mkl.h"


namespace Analysis
{

	void StatusDisplay(long status)
	{
		long classError;                                        
		
		classError = DftiErrorClass(status, DFTI_ERROR_CLASS);  
		if(! classError)
		{
			cerr << "Error Status is not a member of Predefined Error Class\n"; 
		}
		else
		{
			char* errorMessage = DftiErrorMessage(status);
			cerr << "Error_message = " << errorMessage << endl;
		}
	}


	Vector<ComplexFloat> FFT(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}

	Vector<ComplexFloat>& FFTI(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		//Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}




	Vector<ComplexDouble> FFT(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return y;
	}

	Vector<ComplexDouble>& FFTI(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		//Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return x;
	}


	Vector<ComplexFloat> FFT(Vector<float>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_REAL, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}
		
		for(int i=1;i<(n+1)/2;i++)
		{
			y[n-i] = conj(y[i]);
		}

		return y;
	}



	Vector<ComplexDouble> FFT(Vector<double>& x)
	{
		int n = x.Length();
		Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_REAL, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}
		
		for(int i=1;i<(n+1)/2;i++)
		{
			y[n-i] = conj(y[i]);
		}

		return y;
	}








	// ////////////////////
	// Backward transform
	// ////////////////////


	Vector<ComplexFloat> IFFT(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}

	Vector<ComplexFloat>& IFFTI(Vector<ComplexFloat>& x)
	{
		int n = x.Length();
		//Vector<ComplexFloat> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}



	Vector<ComplexDouble> IFFT(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return y;
	}




	Vector<ComplexDouble>& IFFTI(Vector<ComplexDouble>& x)
	{
		int n = x.Length();
		//Vector<ComplexDouble> y(n);
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Backward transform
		status = DftiComputeBackward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}

		return x;
	}
















	// /////////////////
	// 2D Transform
	// /////////////////


	Matrix<ComplexFloat> FFT2(Matrix<ComplexFloat>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		Matrix<ComplexFloat> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}


	Matrix<ComplexFloat>& FFT2I(Matrix<ComplexFloat>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		//Matrix<ComplexFloat> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}




	Matrix<ComplexDouble> FFT2(Matrix<ComplexDouble>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		Matrix<ComplexDouble> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}


	Matrix<ComplexDouble>& FFT2I(Matrix<ComplexDouble>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		//Matrix<ComplexDouble> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeForward( DescHandle, x.Data() );
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}












	Matrix<ComplexFloat> IFFT2(Matrix<ComplexFloat>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		Matrix<ComplexFloat> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}


	Matrix<ComplexFloat>& IFFT2I(Matrix<ComplexFloat>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		//Matrix<ComplexFloat> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeBackward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}




	Matrix<ComplexDouble> IFFT2(Matrix<ComplexDouble>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		Matrix<ComplexDouble> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return y;
	}




	Matrix<ComplexDouble>& IFFT2I(Matrix<ComplexDouble>& x)
	{
		int rows = x.Rows();
		int cols = x.Columns();
		//Matrix<ComplexDouble> y(rows,cols);
		int dims[2] = {cols,rows};
		int strides[3] = {0,rows,1};
		
		DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
		long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Strides
		status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		// Not Inplace
		//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		//if(status != DFTI_NO_ERROR)
		//{
		//	StatusDisplay(status);
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Problem with FFT!");
		//}
		// Backward scale
		status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Commit Dfti descriptor
		status = DftiCommitDescriptor(DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}

		// Compute Forward transform
		status = DftiComputeBackward( DescHandle, x.Data());
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Problem with FFT!");
		}
		status = DftiFreeDescriptor(&DescHandle);
		if(status != DFTI_NO_ERROR)
		{
			StatusDisplay(status);
			Utility::Warning("Problem when trying to free the FFT descriptor handle!");
		}


		return x;
	}












//=========================================================================
// Convolution
//=========================================================================


	Vector<float> Conv(Vector<float>& input, Vector<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
	}


	Vector<double> Conv(Vector<double>& input, Vector<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
	}



	Vector<ComplexFloat> Conv(Vector<ComplexFloat>& input, Vector<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter size is smaller than the input size
		if(filter.Length() > input.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter length should not be larger than the input signal!");
		}

		int fftlength = input.Length()+filter.Length()-1;
		if(PowerOfTwo)
		{
			fftlength = NextPow2(fftlength);
		}
		

		Vector<ComplexFloat> temp1;
		if(b == ZeroPad || b == Symmetric || b == Replicate) 
		{
			temp1 = Vector<ComplexFloat>(fftlength);
			for(int i=0; i<input.Length(); i++)
			{
				temp1[i] = input[i];
			}
		}
		else if(b == Circular)
		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}
		
		
		if(b == Symmetric)
		{
			int wid1 = filter.Length()/2;
			int wid2 = filter.Length() - wid1 - 1;
			if(wid1 != 0)
			{
				Vector<ComplexFloat> c1 = Reverse( input.Slice(0, wid1-1) );
				temp1.ReadFromVector(c1, temp1.Length()-wid1);
			}
			if(wid2 != 0)
			{
				Vector<ComplexFloat> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
				temp1.ReadFromVector(c2, input.Length());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Length()/2;
			int wid2 = filter.Length() - wid1 - 1;
			if(wid1 != 0)
			{
				Vector<ComplexFloat> c1(wid1, input.ElemNC(0));
				temp1.ReadFromVector(c1, temp1.Length()-wid1);
			}
			if(wid2 != 0)
			{
				Vector<ComplexFloat> c2(wid2, input.ElemNC(input.Length()-1));
				temp1.ReadFromVector(c2, input.Length());
			}
		}

		Vector<ComplexFloat> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Vector<ComplexFloat>(fftlength);
			for(int i=0; i<filter.Length(); i++)
			{
				temp2[i] = filter[i];
			}
		}
		else if(b == Circular)
		{
			temp2 = Vector<ComplexFloat>(input.Length());
			for(int i=0; i<filter.Length(); i++)
			{
				temp2[i] = filter[i];
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}

		
		Vector<ComplexFloat> o = IFFT( FFT(temp1) * FFT(temp2) );
		
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Length()+filter.Length()-2);
		}
		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;
		
	}

	Vector<ComplexDouble> Conv(Vector<ComplexDouble>& input, Vector<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter size is smaller than the input size
		if(filter.Length() > input.Length())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter length should not be larger than the input signal!");
		}

		int fftlength = input.Length()+filter.Length()-1;
		if(PowerOfTwo)
		{
			fftlength = NextPow2(fftlength);
		}
		

		Vector<ComplexDouble> temp1;
		if(b == ZeroPad || b == Symmetric || b == Replicate) 
		{
			temp1 = Vector<ComplexDouble>(fftlength);
			for(int i=0; i<input.Length(); i++)
			{
				temp1[i] = input[i];
			}
		}
		else if(b == Circular)
		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}
		
		
		if(b == Symmetric)
		{
			int wid1 = filter.Length()/2;
			int wid2 = filter.Length() - wid1 - 1;
			if(wid1 != 0)
			{
				Vector<ComplexDouble> c1 = Reverse( input.Slice(0, wid1-1) );
				temp1.ReadFromVector(c1, temp1.Length()-wid1);
			}
			if(wid2 != 0)
			{
				Vector<ComplexDouble> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
				temp1.ReadFromVector(c2, input.Length());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Length()/2;
			int wid2 = filter.Length() - wid1 - 1;
			if(wid1 != 0)
			{
				Vector<ComplexDouble> c1(wid1, input.ElemNC(0));
				temp1.ReadFromVector(c1, temp1.Length()-wid1);
			}
			if(wid2 != 0)
			{
				Vector<ComplexDouble> c2(wid2, input.ElemNC(input.Length()-1));
				temp1.ReadFromVector(c2, input.Length());
			}
		}

		Vector<ComplexDouble> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Vector<ComplexDouble>(fftlength);
			for(int i=0; i<filter.Length(); i++)
			{
				temp2[i] = filter[i];
			}
		}
		else if(b == Circular)
		{
			temp2 = Vector<ComplexDouble>(input.Length());
			for(int i=0; i<filter.Length(); i++)
			{
				temp2[i] = filter[i];
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}

		
		Vector<ComplexDouble> o = IFFT( FFT(temp1) * FFT(temp2) );
		
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Length()+filter.Length()-2);
		}
		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;
		
	}
	
	
	// =============
	// Conv2
	// =============


	Matrix<float> Conv2(Matrix<float>& input, Matrix<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv2(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
	}

	Matrix<float> Conv2(Vector<float>& filter1, Vector<float>& filter2, Matrix<float>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv2(ToComplexFloat(filter1), ToComplexFloat(filter2), ToComplexFloat(input), fa, b, PowerOfTwo) );
	}

	Matrix<double> Conv2(Matrix<double>& input, Matrix<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv2(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
	}

	Matrix<double> Conv2(Vector<double>& filter1, Vector<double>& filter2, Matrix<double>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		return Real( Conv2(ToComplexDouble(filter1), ToComplexDouble(filter2), ToComplexDouble(input), fa, b, PowerOfTwo) );
	}


	
	Matrix<ComplexFloat> Conv2(Matrix<ComplexFloat>& input, Matrix<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{

		// Check if the filter size is smaller than the input size
		if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
		}
		
		
		// Allocate area for input using the border
		
		int fftrows = input.Rows()+filter.Rows()-1;
		int fftcols = input.Columns()+filter.Columns()-1;
		if(PowerOfTwo)
		{
			fftrows = NextPow2(fftrows);
			fftcols = NextPow2(fftcols);
		}
		
		Matrix<ComplexFloat> temp1;
		if(b == ZeroPad || b == Symmetric || b == Replicate) 
		{
			temp1 = Matrix<ComplexFloat>(fftrows, fftcols,0);
			for(int j=0; j<input.Columns(); j++)
			{
				for(int i=0; i<input.Rows(); i++)
				{
					temp1.ElemNC(i,j) = input.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		if(b == Symmetric)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexFloat> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
				temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexFloat> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
				temp1.ReadFromMatrix(c2, 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexFloat> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexFloat> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r2, input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexFloat> c1 = input.Slice(0, input.Rows()-1, 0, 0);
				temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexFloat> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexFloat> r1 = input.Slice(0, 0, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexFloat> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q11( hei1, wid1, input.ElemNC(0,0) );
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexFloat> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexFloat> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}

		
		// Allocate area for the filter
		Matrix<ComplexFloat> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Matrix<ComplexFloat>(fftrows, fftcols,0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp2 = Matrix<ComplexFloat>(input.Rows(), input.Columns(),0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		// FFT and inverse FFT
		Matrix<ComplexFloat> o = IFFT2( FFT2(temp1)*FFT2(temp2) );

		// Based on the filter area, crop the border.
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
		}

		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;
		
	}


	
	Matrix<ComplexFloat> Conv2(Vector<ComplexFloat>& filter1, Vector<ComplexFloat>& filter2, Matrix<ComplexFloat>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter sizes is smaller than the input size
		if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
		}
		
		Vector<ComplexFloat> *cols = new (std::nothrow) Vector<ComplexFloat>[input.Columns()];
		Utility::CheckPointer(cols);
		
		for(int i=0; i<input.Columns(); i++)
		{
			cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
		}
		
		int rlength = cols[0].Length();
		
		Vector<ComplexFloat> *rows = new (std::nothrow) Vector<ComplexFloat>[rlength];
		Utility::CheckPointer(rows);
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Vector<ComplexFloat>(input.Columns());
			for(int j=0; j<input.Columns(); j++)
			{
				rows[i][j] = cols[j][i];
			}
		}
		
		delete [] cols;
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
		}
		
		int clength = rows[0].Length();
		Matrix<ComplexFloat> temp(rlength, clength);
		for(int i=0; i<rlength; i++)
		{
			for(int j=0; j<clength; j++)
			{
				temp.ElemNC(i,j) = rows[i][j];
			}
		}
		
		delete [] rows;
		
		return temp;
	}




	Matrix<ComplexDouble> Conv2(Matrix<ComplexDouble>& input, Matrix<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter size is smaller than the input size
		if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
		}
		
		
		// Allocate area for input using the border
		
		int fftrows = input.Rows()+filter.Rows()-1;
		int fftcols = input.Columns()+filter.Columns()-1;
		if(PowerOfTwo)
		{
			fftrows = NextPow2(fftrows);
			fftcols = NextPow2(fftcols);
		}
		
		Matrix<ComplexDouble> temp1;
		if(b == ZeroPad || b == Symmetric || b == Replicate) 
		{
			temp1 = Matrix<ComplexDouble>(fftrows, fftcols,0);
			for(int j=0; j<input.Columns(); j++)
			{
				for(int i=0; i<input.Rows(); i++)
				{
					temp1.ElemNC(i,j) = input.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp1 = input;
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}
		
		if(b == Symmetric)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexDouble> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
				temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexDouble> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
				temp1.ReadFromMatrix(c2, 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexDouble> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexDouble> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
				temp1.ReadFromMatrix(r2, input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}
		else if(b == Replicate)
		{
			int wid1 = filter.Columns()/2;
			int wid2 = filter.Columns() - wid1 - 1;
			if(wid1 != 0)
			{
				Matrix<ComplexDouble> c1 = input.Slice(0, input.Rows()-1, 0, 0);
				temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
			}
			if(wid2 != 0)
			{
				Matrix<ComplexDouble> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
			}
			
			int hei1 = filter.Rows()/2;
			int hei2 = filter.Rows() - hei1 - 1;
			if(hei1 != 0)
			{
				Matrix<ComplexDouble> r1 = input.Slice(0, 0, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
			}
			if(hei2 != 0)
			{
				Matrix<ComplexDouble> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
				temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
			}
			
			
			if(hei1 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q11( hei1, wid1, input.ElemNC(0,0) );
				temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
			}

			if(hei1 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
				temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
			}

			if(hei2 != 0 && wid1 != 0)
			{
				Matrix<ComplexDouble> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
				temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
			}

			if(hei2 != 0 && wid2 != 0)
			{
				Matrix<ComplexDouble> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
				temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
			}
		}

		
		// Allocate area for the filter
		Matrix<ComplexDouble> temp2;
		if(b == ZeroPad || b == Symmetric || b == Replicate)
		{
			temp2 = Matrix<ComplexDouble>(fftrows, fftcols,0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else if(b == Circular)
		{
			temp2 = Matrix<ComplexDouble>(input.Rows(), input.Columns(),0);
			for(int j=0; j<filter.Columns(); j++)
			{
				for(int i=0; i<filter.Rows(); i++)
				{
					temp2.ElemNC(i,j) = filter.ElemNC(i,j);
				}
			}
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Invalid Border type!");
		}


		// FFT and inverse FFT
		Matrix<ComplexDouble> o = IFFT2( FFT2(temp1) * FFT2(temp2) );
		
		// Based on the filter area, crop the border.
		if(fa == Same && b != Circular)
		{
			o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
		}
		else if(fa == Valid)
		{
			o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
		}
		else if(fa == Full && PowerOfTwo)
		{
			o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
		}
		//else if(fa != Full)
		//{
		//	cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		//	Utility::RunTimeError("Invalid FilterArea type!");
		//}
		
		return o;		
	}
	
	Matrix<ComplexDouble> Conv2(Vector<ComplexDouble>& filter1, Vector<ComplexDouble>& filter2, Matrix<ComplexDouble>& input, FilterArea fa, Border b, bool PowerOfTwo)
	{
		// Check if the filter sizes is smaller than the input size
		if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
		}
		
		Vector<ComplexDouble> *cols = new (std::nothrow) Vector<ComplexDouble>[input.Columns()];
		Utility::CheckPointer(cols);
		
		for(int i=0; i<input.Columns(); i++)
		{
			cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
		}
		
		int rlength = cols[0].Length();
		
		Vector<ComplexDouble> *rows = new (std::nothrow) Vector<ComplexDouble>[rlength];
		Utility::CheckPointer(rows);
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Vector<ComplexDouble>(input.Columns());
			for(int j=0; j<input.Columns(); j++)
			{
				rows[i][j] = cols[j][i];
			}
		}
		
		delete [] cols;
		
		for(int i=0; i<rlength; i++)
		{
			rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
		}
		
		int clength = rows[0].Length();
		Matrix<ComplexDouble> temp(rlength, clength);
		for(int i=0; i<rlength; i++)
		{
			for(int j=0; j<clength; j++)
			{
				temp.ElemNC(i,j) = rows[i][j];
			}
		}
		
		delete [] rows;
		
		return temp;
	}





};